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ABSTRACT fJSF 6# ' 

PDF methods provide a means of calculating the properties of turbulent reacting flows. They have been 
successfully applied to many turbulent flames, including some with finite-rate kinetic effects. Here the 
methods are reviewed with an emphasis on computational issues and on their application to turbulent 
combustion. 
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1. INTRODUCTION 

PDF methods have been applied to a variety of turbulent flows both with, and without, combustion. 
Generally, single-phase, low-Mach-number flows have been considered, in which radiation is not a major 
factor. For such flows, the fundamental dependent variables are the velocities U(x,t) and the compositions 

<k(x»t) (e.g., the species mass fractions and enthalpy). Different probabilistic approaches to modelling 

turbulent flows can be categorized according to the statistics of U(x,t) and ^(x,t) that are considered. For 

example, in a mean-flow closure <U(x,t)> and <(£>(x,t)> are the primary dependent variables: in second- 

t « t 

order closures the variances <u;Uj>, <4> a 4 ) p> and covariances <ui<j) a > are also included. Angled brackets 

denote means (i.e., mathematical expectations) and uj(x,t) = Ui(x,t)-<Uj(x,t)> and (J> a (x,t) = <|> a (x,t) 

-<0a(x,t)> are the fluctuating components of U; and (j) a , respectively. 

In pdf methods, the dependent variable is a probability density function (pdf) or joint pdf of U(x,t), 
(j)(x,t). The pdf contains information equivalent to all the moments; and hence pdf methods are, in this sense, 
more comprehensive than moment closures (e.g. second-order closures). The methods that have proved 
most successful are based on one-point, one-time pdfs, which contain information at each point in the flow 
separately, but no joint information at two or more distinct points. 

In the last 15 years, pdf methods have advanced from being only of theoretical interest to a small group 
of specialists, to being a practical approach for calculating the properties of turbulent reactive flows. As the 
review below indicates, in additional to having been applied to idealized flames and simple laboratory flames, 
the methods have been applied to flames requiring multistep chemical kinetics (e.g. Pope 1981a, Chen & 
Kollmann 1988a), as well as to computationally difficult flows (e.g. that in the cylinder of a spark-ignition 
engine, Haworth & El Tahry 1989a, b). 

The purpose of this paper is to review the work on pdf methods, with some emphasis on the numerical 
issues, and on the applications to turbulent combustion. 

In the next section the different pdf methods are described, along with the associated modelling. Monte 
Carlo methods have proved to be the most successful means of solving pdf transport equations. The essence 
of these solution techniques is described in Section 3. The theoretical foundations of pdf methods (including 
the modelling and Monte Carlo solution algorithms) are comprehensively described by Pope (1985). The 
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purpose of Sections 2 and 3 is to describe briefly the principal features, with no attempt at rigor. Section 4 

reviews the applications of pdf methods to turbulent diffusion flames and premixed flames. (Recent 
applications to constant-density inert flows have been reviewed by Pope 1989a.) In the Discussion and 
Conclusions some of the outstanding problems and future directions are assessed. 

2 . PDF METHODS 


2.1 Definitions and Properties 

Let (j) denote the value of a composition variable (the mass fraction of oxygen, for example) at a particular 
location x 0 and time t Q in a turbulent reactive flow. It is supposed that the flow can be realized any number 
of times, and the time t is measured from the initiation of the flow. Thus from each realization we obtain a 

value of §; and, given the nature of turbulence, these values are, in all probability, different. In other words 

(j) is a random variable. On a given realization, it is not possible to determine beforehand the value of <|) that 
will obtain. But it is possible to ascribe probabilities to its value being in a given interval: this can be done 
through the pdf. 

For every random variable we introduce an independent (sample-space) variable: in particular, \j / is the 
sample-space variable corresponding to (j>. The cumulative distribution function (cdf), F<j)(\j/), is then defined 
as the probability that (j) is less than \|/: 


F<|)(\|/) ^ Prob{<j><\|/} . (1) 

And the pdf of <j>, f^tj/), is defined by 

W^F#- (2) 

dvj/ 

While F<j)(\|/) is a probability function, f(j>(\}/) is a probability density function. That is, f<j,(\|/) is the 
probability per unit y/ of the event <| >=\j/ : or, equivalently, f < j ) (\j/)d\)/ is the probability of the event 

¥<<{><¥+d¥- 

The three fundamental properties of the pdf (in addition to Eq. 2) are: 

f<f>(¥) > 0 , (3) 

since probabilities are non-negative; 

oo 

j f^Y) d\|/=l , (4) 

-oo 

since Prob { <{)<oo } = 1 and Prob { c()< — <=<> ) = 0; and, for any (non-pathological) function Q(<|>): 

OO 

<Q> = J f<t>(¥) Q (¥)d¥ • (5) 
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This last equation shows that if the pdf is known, then the expectation of any function of the random variable 
can be calculated. In particular the mean <(j)> and the rn-th central moment <4>' m > (m>l) can be determined 
(if it exists). 

For a general turbulent reactive flow we need to consider a set of a > 1 composition variables (J) = 
{01>02, ...,<j)o}- Accordingly the O' sample-space variables {\}/l>Y2> are introduced, and the 

joint pdf of <[>, f^(\jr), is defined to be the probability density of the compound event <j) = \j/ (i.e. (j>i = \j/i, (j>2 = 
Y2. 0a = Va). 

Clearly the joint pdf defined at the particular location Xo and time t 0 , can be defined at any (x,t): we 
denote by %(xj£;x,t) the joint pdf of <J)(x,t). It is important to realize that this is a one-point, one-time joint 
pdf: it contains no joint information between <]) at two or more positions or times. The pdf method described 
in the next sub-section is based on f^(\|/;x,t), which is called the composition joint pdf. 

The second pdf method described (in Section 2.3) is based on the velocity-composition joint pdf, 
f(V,\jr,x,t). Here, V = {Vi,V 2 ,V 3 } are the three independent velocity variables; and f is the probability 
density of the compound event {U(x,t) = V, 0(x,t) = W)- 

In the treatment of variable-density flows, two further probability functions prove useful, and are now 
defined. By assumption (see Pope 1985), the composition variables are sufficient to determine the fluid 

density. If the composition is 0, the density is given by the function p 0 ((j>), which can be determined from a 
thermodynamic calculation. Consequently at (x,t) the fluid density is 

P(x,t) = pa(0[x>t]) • (6) 

Having made the distinction between the different functions p(x,t) and p o (0), we now follow conventional 
(if imprecise) notation and denote both by p. 

Favre, or density-weighted, pdfs are defined by, for example: 

7<t>(¥) = P(V)f<t>(Y)/<P> • (7) 

It then follows that density-weighted means are given by 

oo 

Q- <pQ>/<p> = j Q(\|/)?(\|0d\|/ , (8) 


(cf. Eq. 5). The mass density function F is defined by 

E(V,y/,x;t) = p(n/)f(V,\j/;x,t) . (9) 

The use of these functions is made apparent in the next two subsections. 


337 



2.2 Composition Joint PDF Equation 


Dopazo and O'Brien (1974) were the first to consider the transport equation for %(\j/;x,t). Since then a 
number of different derivations have been given (e.g. Pope 1976, Janicka et al. 1978a, b, O'Brien 1980, 
Pope 1985). Here we state the result, and refer the reader to Pope (1985) for a detailed derivation. 

The compositions evolve according to the conservation equation 


Dfoq 

Dt 


3J a 

+Sa • 

p 3xi 


( 10 ) 


where J a is the (molecular) diffusive flux of (|) a , and S a — -a known function of (j) — is the rate of creation 
of <f>a due to chemical reaction. The pdf transport equation corresponding to Eq. (10) is 


3t 1 3xi 


+ r~~ [S a (S|£)^ 
3Va 


1 3 

<p> 3xi 


[<uJ\j£>T^] 


3 

+ 

3Va 




— • 
C/Xj 


( 11 ) 


On the left-hand side, the first two terms represent the rate of change following the mean flow. The third 
term is — in composition space — the divergence of the flux of probability due to reaction. It is the form of 

this term that gives this pdf method the advantage over other statistical approaches: since S(\|/) is known, 

is the subject of the equation, and \)/ a is an independent variable, the term contains no unknowns. Thus 
however complicated and non-linear the reaction scheme, in the composition joint pdf equation the effect of 
reaction is in closed form, requiring no modelling. 

In contrast the terms on the right-hand side require modelling. The quantity <u"l\j/> is the mean of the 

Favre velocity fluctuation (u" = U-U) conditional upon the event ^ = rj/. It represents the transport of f^ in 
physical space by the fluctuating velocity. Although there have been other suggestions, this term is generally 
modelled by gradient diffusion: 


<P><Uj I w > *4 = - r T — 

C/Xj 



( 12 ) 


where Fp is a turbulent diffusivity. Such gradient transport models are, of course, subject to many 
objections, especially when applied to variable-density reactive flows. 

The final term in Eq. (11) represents the effect of molecular mixing. It is generally treated by a stochastic 
mixing model (see e.g. Flagan & Appleton 1974, Pope 1982, 1985). While some aspects of this modelling 
are discussed below, the cited references should be consulted for a full account. 

The composition joint pdf equation is not a self-contained model: mean momentum equations must be 
solved for U; and a turbulence model (k-e, say) is needed to determine both Tp (~ k 2 /e) and the mixing rate 
(~ e/k) used in the stochastic mixing model. 
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2.3 Velocity-Composition Joint PDF Equation 

Two shortcomings of the composition pdf approach are that turbulent transport (<u.ljj/>) has to be 

modelled, and that the velocity and turbulence fields have to be treated separately. Both these shortcomings 
are overcome in the velocity-composition joint pdf approach. 

The instantaneous momentum equation is 


DU; 3p 

p ^ = ^‘^ +psi 


where tij is the stress tensor, p the pressure, and g the gravitational acceleration. 


(13) 


From this equation (and that for <j), Eq. 10) the following equation can be derived (Pope 1985) for the 
mass density function F(V,\j/,x;t) (Eq. 9): 

dF T7 dF r _ 1 3<p> 7 dF d rT _, 0 

— + Vj — + [gj - —— — ] + [F S a (ij/)] - 

3t 3xj p(\j/) 3xj 3Vj d\\f a 


3J a 

A.[<.^iL + ^IV,^F/p(v|/)]+-^- [<-- lV,\j£> F/p(\j/>] . 

aVj oxi aXj o\\f(x c/Xj 


(14) 


None of the terms on the left-hand side requires modelling. In order, the terms represent: rate of change 
with time; transport in position space (by both mean and fluctuating components of velocity); transport in 
velocity space (by gravity and the mean pressure gradient); and, as before, transport in composition space 
due to reaction. 


The terms requiring modelling (on the right-hand side of Eq. 1) are means conditional on the compound 

event {U(x,t) = V, 4>(x,t) = \^}. The final term — as in the composition pdf equation — represents 
molecular mixing. The remaining term represents transport in velocity space due to molecular stresses and 
the fluctuating pressure gradient. A discussion of its modelling is deferred to the next subsection. 

It may be seen, then, that the velocity-composition joint pdf method retains the advantage of treating 
reaction without approximation, and, in addition, treats transport in physical space (turbulent convection) 
exactly, thus avoiding gradient-diffusion assumptions. It also provides a more complete closure: The mean 

velocity U(x,t), the Reynolds stresses, and indeed all one-point velocity-composition statistics can be 
calculated from F. The model equation for F is not quite self-contained because the modelled terms require a 

knowledge of the turbulent time scale (k/e) which cannot be deduced from F. 

2.4 Lagrangian Viewpoint 


Thus far the Eulerian view has been adopted: we have considered functions (e.g. F(V,\|/,x;t)) at a fixed 
position x. It proves extremely helpful, both to the modelling and to the numerical solution technique, to 
take the alternative Lagrangian viewpoint also. 


Let x+(t), U + (t) and ^ + (t) denote the position, velocity and composition of the fluid particle that was at a 
reference point Xq at a reference time to. These particle properties evolve according to: 
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dx + (t) 

dt 


U + (t) = U(x + [t],t), 


(15) 


since, by definition, a fluid particle moves with the local fluid velocity; 


du; 

dt 


A i = gj- 


3<p> 
p(^+) L 3xj . 


+ &. ii. -J— ^1 

x+ l9xi p(^ + ) 3xjJx+ 


(from Eq. 13); and, 


d<}> 


a 


"df = e o sS a 0fc + ) 


p(^ + ) 


3J a 

1 


3xJ 


(16) 


(17) 


(from Eq. 10). 

The connection between these equations for the properties of a fluid particle, and the equation for the 
mass density function F (Eq. 14) is immediately apparent. Equation (14) can be written 


— + — [F<U + I>] + — [F<A;I>] + — [F<0 a l>] = 0, (18) 

3t 3xj J 3Vj 3\|/a 


where the expectations are conditional on the compound event (x + (t) = x, U + (t) = V, (j> + (t) = tj/}. Further, it 
may be noticed that the terms in braces in Eqs. (16)-(17) appear on the right-hand side of Eq. (14) — that is, 
they need to be modelled — whereas all other terms appear on the left-hand side and are treated exactly. 

2.5 Stochastic Models 

The standard approach to turbulence modelling is to construct constitutive relations for the unknown 
correlations (see, e.g. Lumley 1978). In the context of the mass density function, this approach is to model 
the unknown conditional expectations on the right-hand side of Eq. (14) in terms of known quantities, i.e. 

functions or functionals of F(V,\j/,x;t). But the Lagrangian viewpoint offers a different approach to 

modelling: namely to use stochastic processes to simulate unknown contributions to U+(t) and (j) + (t) (i.e. the 
terms in braces in Eqs. 16 and 17). 

To illustrate this approach we consider U*(t) — a stochastic model for U + (t). If the model is accurate 
then U*(t) is (statistically) an accurate approximation to U + (t). In general the time series U* is not 
differentiable. Consequently we express the models in terms of the infinitesimal increment 


dU*(t) = n*(t+dt) - U*(t) , (19) 

rather than in terms of the derivative dU*/dt. Note that for a deterministic, differentiable process, (e.g. 
U + (t)), the infinitesimal increment is non-random (i.e. zero variance) and is of order dt. 

In view of the equation for U + (t) (Eq. 16), the increment dU* can be written 



f 1 

d<p> 


l 6j ’ pm 

dxj 

x* y 


dt + dUj , 


( 20 ) 
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where (similar to U*) x* and <|)* are models of x + and (J> + . The stochastic increment dU s models the effects 
of the fluctuating pressure gradient and viscous stresses, while the term in dt is an exact expression for the 
effect of gravity and the mean pressure gradient. 

Two types of models have been used. The first type — of which the stochastic mixing model is an 
example — are called particle interaction models. In the terminology of stochastic processes, these are point 
processes. According to these models, the infinitesimal increment dU s is nearly always zero. But with 
probability of order dt, the increment is of order unity. Thus the time series is a piecewise constant, with of 
order unity jumps per unit time. 

The second type of model use diffusion processes in which dU s is a random variable with (conditional) 
mean and variance both of order dt. Note that this implies that the rms is of order dt 1 / 2 , and hence the 
process — though continuous — is not differentiable. The different variants of the Langevin model are 
diffusion processes (see e.g. Pope 1983, 1985, Haworth & Pope 1986, 1987a). 

For more information on this general modelling approach the reader is referred to Pope (1985), while the 
current status of the Langevin model is described by Haworth & Pope (1987a). 

3 . NUMERICAL SOLUTION ALGORITHMS 

The velocity-composition joint pdf f(V,\j/;x,t) is a single function defined in a multi-dimensional space. 
In general f depends on the three velocity variables, a composition variables, three spatial variables, and time 

— (7+g) independent variables in all. In many cases the dimensionality may be less, but still large. For 
example, in a statistically stationary and two-dimensional flow with a single composition variable, 
f(Y>YF x l> x 2) depends on 6 independent variables. The composition joint pdf f^(\j/;x,t) in general depends 
on (4+a) variables; but, for the simpler flow cited above, f<f>(>|/i;xi,x2) is a function of just 3 variables. 

Given the large dimensionality of joint pdfs it is clear that conventional grid-based numerical methods 
(e.g. finite differences) are impractical for all but the simplest of cases. Just to provide an accurate 
representation of a function of 6 independent variables is a major task. Consequently, although one or two 

finite-difference solutions have been obtained for f<j>(vj/i;xi,x 2 ) (e.g. Janicka & Kollmann 1978a,b), currently 
all investigators are using Monte Carlo methods instead. 

In the next subsection the general Monte Carlo method devised by Pope (1985) to solve for the velocity- 
composition joint pdf is outlined. Then, in section 3.2, Monte Carlo solution algorithms for the composition 
joint pdf are reviewed. 

3.1 Monte Carlo Method for the Velocity-Composition Joint PDF 

The Monte Carlo method to solve the modelled equation for the velocity-composition joint pdf is 
conceptually simple and natural. Rather than discretizing the space, we discretize the mass of fluid into a 
large number N of representative or stochastic particles. At a given time t, let M be the total mass of fluid 

within the solution domain. Then each stochastic particle represents a mass Am = M/N of fluid. The n-th 

particle has position x( n )(t), velocity lK n )(t) ; and composition <j>( n )(t). 

Starting from appropriate initial conditions, the particle properties are advanced in time by the increments 

dx( n )(t) = U( n )(t) dt , (21) 

dU(“)(t) = [g-p(^ n ))' 1 V<p>]dt + dU s , (22) 
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d(j)( n )(t) = S(§( n ))dt + d(j) s , 


(23) 


where dU s and d<j> s are the stochastic increments that simulate molecular processes and the fluctuating 
pressure gradient. At symmetry boundaries particles are reflected; at inflow boundaries particles are added 
with appropriate properties; and, at outflow boundaries particles are discarded. While wall boundaries have 
been treated (Anand et al. 1989a) a comprehensive account of this treatment is not available in the literature. 


The correspondence between the ensemble of stochastic particles and the joint pdf has been established 
by Pope (1985). The main results are: 

N 

i) The expected density of the stochastic particles in physical space (Am £ <8(x-x( n )(t)>) is equal to the 

n=l 


fluid density <p(x,t)>. 


ii) 


The joint pdf of the stochastic particle properties U( n )(t), <|)( n )(t) is the density-weighted joint pdf 

7(V,\|/;x(n)[t],t). 


iii) From particle properties, expectations (e.g. U(x,t)) can be approximated as ensemble averages, with a 
statistical error of order N' 1 / 2 . 


Several implementations of this algorithm, and variants of it, have been reported. For example, the 
turbulent jet diffusion flame calculations reported in Section 4 are performed using a "boundary-layer" 
variant (see Pope 1985). Haworth & Pope (1987b) report a variant of the algorithm designed specifically for 
self-similar shear flows. From a numerical standpoint this work is of particular interest, because the 
convergence of the method (as N' 1 / 2 ) is demonstrated. The basic algorithm has been implemented and 
demonstrated for statistically two-dimensional recirculating flows by Anand et al. (1989b). Haworth & El 
Tahry (1989a,b) reports calculations of the three-dimensional time-dependent flow in the cylinder of a spark- 
ignition engine. In this case the pdf algorithm is coupled to a conventional finite-volume algorithm that is 
used to calculate the mean pressure and turbulent time scale fields. 

3.2 Monte Carlo Algorithms for the Composition Joint PDF 

Two different algorithms have been proposed to solve the modelled transport equation for the 
composition joint pdf. 

The second of these (chronologically) was proposed by Pope (1985), and is similar to that described 
above for the velocity-composition joint pdf. Again, it is a grid-free algorithm in which the mass of fluid is 

discretized into N stochastic particles, the n-th of these having position x( n )(t) and composition ^( n )(t). In 
each time step the composition is incremented according to Eq. (23), while the position is incremented by 


dx( n )(t) = D(x (n) [t],t)dt + dx s , (24) 

where the stochastic component dx s causes a random walk to simulate gradient diffusion. No 
implementations of this algorithm have been reported in the literature. 

The first Monte Carlo algorithm for the composition joint pdf was devised by Pope (1981b). In this 
method there is a finite-difference grid in physical space. At each grid node, the composition joint pdf is 

represented by N particles, the n-th having composition £( n )(t). Reaction and mixing are performed 
according to Eq. (23), while particles are moved from node to node to simulate convection and turbulent 
diffusion. 
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This algorithm is used in the premixed flame calculation of Pope (1981a) and McNutt (1981), and in the 

diffusion flame calculations of Nguyen & Pope (1984), Jones & Kollmann (1987) and Chen & Kollmann 
(1988a). 


4. TURBULENT FLAME CALCULATIONS 

4.1 Turbulent Diffusion Flames 

Some of the first pdf calculations are of turbulent diffusion flames (Frost 1975, Janicka, Kolbe and 
Kollmann 1978a, b, Bywater 1982, Nguyen & Pope 1984). The calculations reported by Nguyen & Pope 
(1984) are the first use of the Monte Carlo method for jet flames. The results include demonstrations of 
convergence of the solutions as N‘ 1</2 tends to zero. 

In the calculations cited above, the thermochemistry is handled in a simple manner — by assuming 
chemical equilibrium, for example. This reduces the number of compositions to just one; namely, the 
mixture fraction. Finite-rate, multi-step kinetics have been used by Pope & Correa (1986) (see also Correa, 
Gulati & Pope 1988), Jones & Kollmann (1987) and Chen & Kollmann (1988a, b). A computational 
challenge is to implement the integration of the rate equation, i.e. 


-£ = S(<J>) . (25) 

in an efficient manner. All investigators have used table-look-up algorithms. 

Considerable attention has been paid to the CO/H 2 -air turbulent diffusion flame studied experimentally by 
Drake et al. (1984). Using the velocity-composition joint pdf approach Pope & Correa (1986) and Correa, 
Gulati and Pope (1988) report calculations based on a partial equilibrium assumption. This reduces the 
number of compositions to two: the mixture fraction and a reaction progress variable (for the radical 
recombination reactions). Again using the velocity-composition joint pdf approach, Haworth, Drake & Blint 
(1988) and Haworth, Drake, Pope & Blint (1988) have made calculations of this flame using a flamelet 
model. 

4.2 Turbulent Premixed Flames 

The composition joint pdf approach (using the Monte Carlo method) has been applied to premixed flames 
by Pope (1981a) and McNutt (1981). The purpose of the former calculation was to demonstrate the ability 
of the pdf method to handle nonlinear reaction kinetics. A three-variable kinetic scheme was used to calculate 
the oxidation of CO and formation of NO in a propane-air flame stabilized behind a perforated plate. 

The works of McNutt (1981), Pope & Anand (1984) and Anand & Pope (1987) are concerned with the 
idealized case of a statistically steady and one-dimensional turbulent premixed flame. In the last of these, the 
velocity-composition joint pdf method is used, and the effects of variable density are studied. It is shown 
that, like the Bray-Moss-Libby model (Bray, Libby & Moss 1985), the pdf method is capable of accounting 
for counter-gradient transport and large turbulence energy production due to heat release. The application of 
the method to a spark-ignited turbulent flame ball is described by Pope & Cheng (1986). 

Turbulent premixed combustion usually occurs in the flamelet regime (Pope 1987). This fact presents a 
challenge to any statistical approach, since the small scales of the composition fields are no longer governed 
by the turbulent straining motions, but are determined by reaction and diffusion occurring in thin flame 
sheets. Pope & Anand (1984) present and demonstrate a model applicable to the flamelet regime. But, as 
discussed by Pope (1985, 1987), this model is not entirely satisfactory. 

An alternative approach to treating flamelet combustion is the stochastic flamelet model of Pope & Cheng 
(1988). This can be viewed as a pdf approach, in which a modelled pdf equation is solved by a Monte Carlo 
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method. In this case, however, the pdf is not that of fluid properties (i.e. velocity and composition) but is 
rather the pdf of flamelet properties (i.e. position, area and orientation of flamelets). 

5. DISCUSSION AND CONCLUSION 

The works reviewed in the previous section demonstrate that pdf methods provide a practicable means of 
calculating the properties of turbulent reactive flows. Calculations have been made with thermochemical 
schemes involving up to three composition variables with finite-rate kinetics (e.g. Pope 1981a, Chen & 
Kollmann 1988a). And the Monte Carlo method used to solve the pdf equations has been implemented for a 
variety of flows including two-dimensional recirculating flows (Anand et al. 1989a) and the three- 
dimensional transient flow in a spark- ignition engine (Haworth & El Tahry 1989a,b). 

The most advanced method considered here is the velocity-composition joint pdf approach. Compared to 
moment closures, this approach has the advantage that reaction can be treated exactly without approximation. 
Compared to the composition joint pdf approach it has the advantages that turbulent transport is treated 
exactly, and that a separate turbulence model is not needed to determine the Reynolds stresses. 

A short-coming of the velocity-composition joint pdf approach is that it does not provide a completely 
self-contained model, in that the turbulence frequency <co> = <e>/k must be determined by separate means. 

For example, in some calculations of simple free shear layers, it has been assumed that <co> is constant 
across the flow, and scales with the mean-flow velocity and length scales (e.g. Haworth & Pope 1987a, 
Pope & Correa 1986). In more complex flows, another approach is to solve the standard model equation for 

<e> (e.g. Haworth & El Tahry 1989a) or, similarly, to solve a modelled equation for <co> deduced from 
those for k and <£> (Anand et al. 1989b). 

A natural extension is to consider f(Y.4/,£;x,t) — the joint pdf of velocity, composition and dissipation. 
This is the probability density function of the compound event (U(x,t) = V, (j)(x,t) = \jf, e(x,t) = Q, where 
e(x,t) is the instantaneous mechanical dissipation. Following some preliminary investigations (Pope & 

Haworth 1986, Pope & Chen 1987) a satisfactory model equation for f(V,\|/,£;x,t) has been developed (Pope 
1989b). The incorporation of dissipation within the pdf allows more realistic and accurate modelling. But 

more important, the single equation for f(V,\j/,^;x,t) provides a completely self-contained model for turbulent 
reactive flows. 

We now discuss three areas in which progress can be expected in the next five years. 

The first area is in the turbulent mixing models. As discussed in Section 2, the stochastic mixing models 
used lead to the composition time series being discontinuous — clearly contrary to the physics of the 
problem. Nevertheless, in spite of their lack of physical appeal, the stochastic models have many advantages 
over alternative suggestions, and, for inert mixing their performance is generally acceptable. But for reactive 
flows, especially in the flamelet regime, their performance is highly suspect. We expect that stochastic 
models will be improved and refined to account better for the microstructure of the composition fields, and 
also to allow mixing and reaction to proceed simultaneously at finite rates. 

The second area of expected progress is in the computational implementation of complex kinetics. When 
the Monte Carlo method is used to solve the joint pdf equation for an inert flow involving o compositions, 
the computer time and storage increases at most linearly with o. But with reaction, in a naive 

implementation, for each particle, on each time step, it is necessary to solve the coupled set of a ordinary 
differential equations 
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dt 


S a (<t>i><h,-,<t>a) , a = l,2,...,o . 


(26) 


The right-hand side (which is a combination of reaction rates) is computationally expensive to evaluate, and, 
as is well known, the set of equations is likely to be stiff. Hence such a naive implementation is 

impracticable for all but the lowest values of o. 

As mentioned in Section 4.1, the alternative approach followed by all investigators is to implement Eq. 
(26) more efficiently through a different table look-up scheme (e.g. Pope & Correa 1986, Jones & Kollmann 
1987). To date this has been done on an ad hoc basis. It is expected that progress will be made in the 
development of a general methodology. 

The third area of expected progress is in the determination of the mean pressure field <p(x,t)> within the 
Monte Carlo algorithm. For thin shear flows, the mean pressure is readily determined by invoking the 
boundary-layer approximations. For statistically-stationary, constant-density, two-dimensional recirculating 
flow, an algorithm to determine <p> has been developed and demonstrated (Anand et al. 1989a). But for the 
general case a computationally efficient and robust algorithm needs to be developed. (In the three- 
dimensional transient calculations of Haworth & El Tahry 1989a, b, the Monte Carlo method is coupled to a 
finite-volume code that determines <p>.) 
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